CAPTAIN workflow continental vs highseas

Author

Théophile L. Mouton

Published

February 12, 2025

Visualise CAPTAIN results

R libraries

Code
library(readr)
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(dplyr)
library(gridExtra)
library(biscale)
library(colorspace)
library(grid)
library(jsonlite)
library(here)

FUSE : conservation prioritiy maps

Budget: 0.3

Code
# Read both RDS files from the Data folder
continental_data_FUSE_03 <- readRDS(here::here("Data/FUSE_continental_full_results_averaged_budget0.3_replicates10.rds"))
high_seas_data_FUSE_03 <- readRDS(here::here("Data/FUSE_full_highseas_results_averaged_budget0.3_replicates10.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform both datasets to sf objects and project
continental_sf_FUSE_03 <- st_as_sf(continental_data_FUSE_03, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

high_seas_sf_FUSE_03 <- st_as_sf(high_seas_data_FUSE_03, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Combine the datasets
combined_sf_FUSE_03 <- rbind(
  mutate(continental_sf_FUSE_03, Region = "Continental Waters"),
  mutate(high_seas_sf_FUSE_03, Region = "High Seas")
)

# Project the world map
world_projected_FUSE_03 <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

# 1. Continental Waters Plot
continental_plot_FUSE_03 <- ggplot() +
  geom_sf(data = continental_sf_FUSE_03, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_FUSE_03, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Conservation Priorities in Continental Waters",
       subtitle = "Index: FUSE, Budget: 0.3, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# 2. High Seas Plot
high_seas_plot_FUSE_03 <- ggplot() +
  geom_sf(data = high_seas_sf_FUSE_03, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_FUSE_03, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Conservation Priorities in High Seas",
       subtitle = "Index: FUSE, Budget: 0.3, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# Combined Plot (modified)
combined_plot_FUSE_03 <- ggplot() +
  geom_sf(data = combined_sf_FUSE_03, 
          aes(color = Priority), 
          size = 0.5, 
          alpha = 0.7) +
  geom_sf(data = world_projected_FUSE_03, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Combined Conservation Priorities",
       subtitle = "Continental Waters and High Seas\nIndex: FUSE, Budget: 0.3, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# Display all three plots
#library(patchwork)
continental_plot_FUSE_03 

Code
high_seas_plot_FUSE_03 

Code
combined_plot_FUSE_03

Budget: 0.1

Code
# Read both RDS files from the Data folder
continental_data_FUSE_01 <- readRDS(here::here("Data/FUSE_full_results_continental_averaged_budget0.1_replicates10.rds"))

high_seas_data_FUSE_01 <- readRDS(here::here("Data/FUSE_full_highseas_results_averaged_budget0.1_replicates10.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform both datasets to sf objects and project
continental_sf_FUSE_01 <- st_as_sf(continental_data_FUSE_01, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

high_seas_sf_FUSE_01 <- st_as_sf(high_seas_data_FUSE_01, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Combine the datasets
combined_sf_FUSE_01 <- rbind(
  mutate(continental_sf_FUSE_01, Region = "Continental Waters"),
  mutate(high_seas_sf_FUSE_01, Region = "High Seas")
)

# Project the world map
world_projected_FUSE_01 <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

# 1. Continental Waters Plot
continental_plot_FUSE_01 <- ggplot() +
  geom_sf(data = continental_sf_FUSE_01, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_FUSE_01, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Conservation Priorities in Continental Waters",
       subtitle = "Index: FUSE, Budget: 0.1, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# 2. High Seas Plot
high_seas_plot_FUSE_01 <- ggplot() +
  geom_sf(data = high_seas_sf_FUSE_01, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_FUSE_01, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Conservation Priorities in High Seas",
       subtitle = "Index: FUSE, Budget: 0.1, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# Combined Plot (modified)
combined_plot_FUSE_01 <- ggplot() +
  geom_sf(data = combined_sf_FUSE_01, 
          aes(color = Priority), 
          size = 0.5, 
          alpha = 0.7) +
  geom_sf(data = world_projected_FUSE_01, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Combined Conservation Priorities",
       subtitle = "Continental Waters and High Seas\nIndex: FUSE, Budget: 0.1, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# Display all three plots
#library(patchwork)
continental_plot_FUSE_01 

Code
high_seas_plot_FUSE_01 

Code
combined_plot_FUSE_01

EDGE2 : conservation prioritiy maps

Budget: 0.3

Code
# Read both RDS files from the Data folder
continental_data_EDGE2_03 <- readRDS(here::here("Data/EDGE2_full_results_continental_averaged_budget0.3_replicates10.rds"))

high_seas_data_EDGE2_03 <- readRDS(here::here("Data/EDGE2_full_highseas_results_averaged_budget0.3_replicates10.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform both datasets to sf objects and project
continental_sf_EDGE2_03 <- st_as_sf(continental_data_EDGE2_03, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

high_seas_sf_EDGE2_03 <- st_as_sf(high_seas_data_EDGE2_03, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Combine the datasets
combined_sf_EDGE2_03 <- rbind(
  mutate(continental_sf_EDGE2_03, Region = "Continental Waters"),
  mutate(high_seas_sf_EDGE2_03, Region = "High Seas")
)

# Project the world map
world_projected_EDGE2_03 <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

# 1. Continental Waters Plot
continental_plot_EDGE2_03 <- ggplot() +
  geom_sf(data = continental_sf_EDGE2_03, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_EDGE2_03, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Conservation Priorities in Continental Waters",
       subtitle = "Index: EDGE2, Budget: 0.3, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# 2. High Seas Plot
high_seas_plot_EDGE2_03 <- ggplot() +
  geom_sf(data = high_seas_sf_EDGE2_03, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_EDGE2_03, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Conservation Priorities in High Seas",
       subtitle = "Index: EDGE2, Budget: 0.3, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# Combined Plot (modified)
combined_plot_EDGE2_03 <- ggplot() +
  geom_sf(data = combined_sf_EDGE2_03, 
          aes(color = Priority), 
          size = 0.5, 
          alpha = 0.7) +
  geom_sf(data = world_projected_EDGE2_03, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Combined Conservation Priorities",
       subtitle = "Continental Waters and High Seas\nIndex: EDGE2, Budget: 0.3, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# Display all three plots
#library(patchwork)
continental_plot_EDGE2_03 

Code
high_seas_plot_EDGE2_03 

Code
combined_plot_EDGE2_03

Budget: 0.1

Code
# Read both RDS files from the Data folder
continental_data_EDGE2_01 <- readRDS(here::here("Data/EDGE2_full_results_continental_averaged_budget0.1_replicates10.rds"))

high_seas_data_EDGE2_01 <- readRDS(here::here("Data/EDGE2_full_highseas_results_averaged_budget0.1_replicates10.rds"))

# Get world map data
world <- ne_countries(scale = "medium", returnclass = "sf")

# Define the McBryde-Thomas 2 projection
mcbryde_thomas_2 <- "+proj=mbt_s"

# Transform both datasets to sf objects and project
continental_sf_EDGE2_01 <- st_as_sf(continental_data_EDGE2_01, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

high_seas_sf_EDGE2_01 <- st_as_sf(high_seas_data_EDGE2_01, coords = c("Longitude", "Latitude"), crs = 4326) %>%
  st_transform(crs = mcbryde_thomas_2)

# Combine the datasets
combined_sf_EDGE2_01 <- rbind(
  mutate(continental_sf_EDGE2_01, Region = "Continental Waters"),
  mutate(high_seas_sf_EDGE2_01, Region = "High Seas")
)

# Project the world map
world_projected_EDGE2_01 <- st_transform(world, crs = mcbryde_thomas_2)

# Create the globe bounding box
globe_bbox <- rbind(c(-180, -90), c(-180, 90), 
                    c(180, 90), c(180, -90), c(-180, -90))

# Create the globe border
globe_border <- st_polygon(list(globe_bbox)) %>%
  st_sfc(crs = 4326) %>%
  st_sf(data.frame(rgn = 'globe', geom = .)) %>%
  smoothr::densify(max_distance = 0.5) %>%
  st_transform(crs = mcbryde_thomas_2)

# Create base theme
my_theme <- theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal",
    legend.box = "vertical",
    legend.margin = margin(t = 20, r = 0, b = 0, l = 0),
    legend.title = element_text(margin = margin(b = 10)),
    panel.background = element_rect(fill = "white", color = NA),
    plot.background = element_rect(fill = "white", color = NA),
    panel.grid = element_blank()
  )

# 1. Continental Waters Plot
continental_plot_EDGE2_01 <- ggplot() +
  geom_sf(data = continental_sf_EDGE2_01, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_EDGE2_01, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Conservation Priorities in Continental Waters",
       subtitle = "Index: EDGE2, Budget: 0.1, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# 2. High Seas Plot
high_seas_plot_EDGE2_01 <- ggplot() +
  geom_sf(data = high_seas_sf_EDGE2_01, aes(color = Priority), size = 0.5, alpha = 0.7) +
  geom_sf(data = world_projected_EDGE2_01, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Conservation Priorities in High Seas",
       subtitle = "Index: EDGE2, Budget: 0.1, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# Combined Plot (modified)
combined_plot_EDGE2_01 <- ggplot() +
  geom_sf(data = combined_sf_EDGE2_01, 
          aes(color = Priority), 
          size = 0.5, 
          alpha = 0.7) +
  geom_sf(data = world_projected_EDGE2_01, fill = "lightgrey", color = "lightgrey", size = 0.1) +
  geom_sf(data = globe_border, fill = NA, color = "black", size = 0.5) +
  scale_color_gradientn(
    colors = c("white", "yellow", "darkblue"),
    values = c(0, 0.5, 1),
    name = "Priority",
    guide = guide_colorbar(barwidth = 20, barheight = 0.5, 
                         title.position = "top", title.hjust = 0.5)
  ) +
  labs(title = "Combined Conservation Priorities",
       subtitle = "Continental Waters and High Seas\nIndex: EDGE2, Budget: 0.1, Replicates: 10",
       x = NULL, y = NULL) +
  my_theme

# Display all three plots
#library(patchwork)
continental_plot_EDGE2_01 

Code
high_seas_plot_EDGE2_01 

Code
combined_plot_EDGE2_01

Congruence tests

Budget: 0.1

Code
# Spearman's rank correlation 
# Calculate correlation between FUSE and EDGE2 priorities
congruence_FUSE_EDGE2_01 <- cor.test(combined_sf_FUSE_01$Priority, 
                                    combined_sf_EDGE2_01$Priority, 
                                    method = "spearman")

# Create a scatterplot to visualize the relationship
congruence_plot_01 <- ggplot(data.frame(
  FUSE = combined_sf_FUSE_01$Priority,
  EDGE2 = combined_sf_EDGE2_01$Priority
), aes(x = FUSE, y = EDGE2)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", color = "red") +
  labs(
    title = "Congruence between FUSE and EDGE2 Priorities",
    subtitle = paste("Spearman's rho =", 
                    round(congruence_FUSE_EDGE2_01$estimate, 3),
                    "\np-value <", 
                    format.pval(congruence_FUSE_EDGE2_01$p.value, digits = 3)),
    x = "FUSE Priority (Budget = 0.1)",
    y = "EDGE2 Priority (Budget = 0.1)"
  ) +
  theme_minimal()

# Display results
print(congruence_FUSE_EDGE2_01)

    Spearman's rank correlation rho

data:  combined_sf_FUSE_01$Priority and combined_sf_EDGE2_01$Priority
S = 1.3942e+14, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.6306786 
Code
congruence_plot_01

Code
# Ovelaps test 
# Create binary maps of high priority areas (≥ 0.8)
high_priority_FUSE_01 <- combined_sf_FUSE_01$Priority >= 0.8
high_priority_EDGE2_01 <- combined_sf_EDGE2_01$Priority >= 0.8

# Calculate overlaps
total_high_FUSE <- sum(high_priority_FUSE_01)
total_high_EDGE2 <- sum(high_priority_EDGE2_01)
overlap_areas <- sum(high_priority_FUSE_01 & high_priority_EDGE2_01)

# Calculate percentages
percent_overlap_FUSE <- (overlap_areas / total_high_FUSE) * 100
percent_overlap_EDGE2 <- (overlap_areas / total_high_EDGE2) * 100

# Print results
cat("Number of high-priority cells:\n",
    "FUSE:", total_high_FUSE, "\n",
    "EDGE2:", total_high_EDGE2, "\n",
    "Overlap:", overlap_areas, "\n\n",
    "Percentage of FUSE high-priority areas that overlap with EDGE2:", round(percent_overlap_FUSE, 1), "%\n",
    "Percentage of EDGE2 high-priority areas that overlap with FUSE:", round(percent_overlap_EDGE2, 1), "%\n")
Number of high-priority cells:
 FUSE: 2130 
 EDGE2: 1951 
 Overlap: 1799 

 Percentage of FUSE high-priority areas that overlap with EDGE2: 84.5 %
 Percentage of EDGE2 high-priority areas that overlap with FUSE: 92.2 %
Code
# Function to calculate hotspot overlaps with significance testing using absolute threshold
calculate_priority_overlap <- function(index1, index2, threshold = 0.8, n_permutations = 999) {
  # Identify high priority areas (>= threshold)
  priority1 <- index1$Priority >= threshold
  priority2 <- index2$Priority >= threshold
  
  # Calculate observed overlap
  Ni <- sum(priority1)  # Number of high priority areas in index1
  Nj <- sum(priority2)  # Number of high priority areas in index2
  NT <- length(priority1)  # Total number of cells
  
  Oo <- sum(priority1 & priority2)  # Observed overlap
  Oe <- (Ni * Nj) / NT  # Expected overlap under independence
  
  # Randomization procedure
  random_overlaps <- numeric(n_permutations)
  for(i in 1:n_permutations) {
    random_priority2 <- sample(priority2)  # Randomly permute second index
    random_overlaps[i] <- sum(priority1 & random_priority2)
  }
  
  # Calculate p-value
  if(Oo > Oe) {
    p_value <- sum(random_overlaps >= Oo) / n_permutations
  } else {
    p_value <- sum(random_overlaps <= Oo) / n_permutations
  }
  
  # Calculate percentage of observed overlap
  percent_overlap <- (Oo / min(Ni, Nj)) * 100
  
  return(list(
    observed_overlap = Oo,
    expected_overlap = Oe,
    total_priority_index1 = Ni,
    total_priority_index2 = Nj,
    percent_overlap = percent_overlap,
    p_value = p_value
  ))
}

# Run the analysis for 0.1 budget level with 0.8 threshold
results_01 <- calculate_priority_overlap(
  combined_sf_FUSE_01,
  combined_sf_EDGE2_01,
  threshold = 0.9
)

# Print results
cat("Results for 0.1 budget level (Priority >= 0.8):\n",
    "\nNumber of high priority areas:",
    "\nFUSE:", results_01$total_priority_index1,
    "\nEDGE2:", results_01$total_priority_index2,
    "\n\nOverlap analysis:",
    "\nObserved overlap:", results_01$observed_overlap,
    "\nExpected overlap:", round(results_01$expected_overlap, 2),
    "\nPercentage overlap:", round(results_01$percent_overlap, 1), "%",
    "\np-value:", format.pval(results_01$p_value, digits = 3))
Results for 0.1 budget level (Priority >= 0.8):
 
Number of high priority areas: 
FUSE: 1812 
EDGE2: 1718 

Overlap analysis: 
Observed overlap: 1632 
Expected overlap: 23.7 
Percentage overlap: 95 % 
p-value: <2e-16

Budget: 0.3

Code
# Calculate correlation between FUSE and EDGE2 priorities at 0.3 budget
congruence_FUSE_EDGE2_03 <- cor.test(combined_sf_FUSE_03$Priority, 
                                    combined_sf_EDGE2_03$Priority, 
                                    method = "spearman")

# Create a scatterplot to visualize the relationship
congruence_plot_03 <- ggplot(data.frame(
  FUSE = combined_sf_FUSE_03$Priority,
  EDGE2 = combined_sf_EDGE2_03$Priority
), aes(x = FUSE, y = EDGE2)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", color = "red") +
  labs(
    title = "Congruence between FUSE and EDGE2 Priorities",
    subtitle = paste("Spearman's rho =", 
                    round(congruence_FUSE_EDGE2_03$estimate, 3),
                    "\np-value <", 
                    format.pval(congruence_FUSE_EDGE2_03$p.value, digits = 3)),
    x = "FUSE Priority (Budget = 0.3)",
    y = "EDGE2 Priority (Budget = 0.3)"
  ) +
  theme_minimal()

# Display results
print(congruence_FUSE_EDGE2_03)

    Spearman's rank correlation rho

data:  combined_sf_FUSE_03$Priority and combined_sf_EDGE2_03$Priority
S = 7.8895e+13, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.7910035 
Code
congruence_plot_03

Code
# Calculate priority overlap 
results_03 <- calculate_priority_overlap(
  combined_sf_FUSE_03,
  combined_sf_EDGE2_03,
  threshold = 0.9
)
# Print results
cat("Results for 0.3 budget level (Priority >= 0.8):\n",
    "\nNumber of high priority areas:",
    "\nFUSE:", results_03$total_priority_index1,
    "\nEDGE2:", results_03$total_priority_index2,
    "\n\nOverlap analysis:",
    "\nObserved overlap:", results_03$observed_overlap,
    "\nExpected overlap:", round(results_03$expected_overlap, 2),
    "\nPercentage overlap:", round(results_03$percent_overlap, 1), "%",
    "\np-value:", format.pval(results_03$p_value, digits = 3))
Results for 0.3 budget level (Priority >= 0.8):
 
Number of high priority areas: 
FUSE: 3470 
EDGE2: 3543 

Overlap analysis: 
Observed overlap: 2709 
Expected overlap: 93.62 
Percentage overlap: 78.1 % 
p-value: <2e-16

FUSE and EDGE2 : conservation prioritiy bivariate maps

Budget: 0.3

Code
# Read all RDS files
edge_continental <- readRDS(here::here("Data/EDGE2_full_results_continental_averaged_budget0.3_replicates10.rds"))
edge_highseas <- readRDS(here::here("Data/EDGE2_full_highseas_results_averaged_budget0.3_replicates10.rds"))
fuse_continental <- readRDS(here::here("Data/FUSE_continental_full_results_averaged_budget0.3_replicates10.rds"))
fuse_highseas <- readRDS(here::here("Data/FUSE_full_highseas_results_averaged_budget0.3_replicates10.rds"))

# Get world map data and set projection
world <- ne_countries(scale = "medium", returnclass = "sf")
mcbryde_thomas_2 <- "+proj=mbt_s"

# Function to process and combine data
process_data <- function(edge_data, fuse_data) {
  combined_data <- edge_data %>%
    rename(EDGE_Priority = Priority) %>%
    left_join(fuse_data %>% rename(FUSE_Priority = Priority),
              by = c("Longitude", "Latitude"))
  
  # Normalize priorities to 0-1 range
  combined_data <- combined_data %>%
    mutate(
      EDGE_Priority_Norm = (EDGE_Priority - min(EDGE_Priority)) / (max(EDGE_Priority) - min(EDGE_Priority)),
      FUSE_Priority_Norm = (FUSE_Priority - min(FUSE_Priority)) / (max(FUSE_Priority) - min(FUSE_Priority))
    )
  
  # Transform to sf object
  data_sf <- st_as_sf(combined_data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
    st_transform(crs = mcbryde_thomas_2)
  
  return(data_sf)
}

# Process continental and high seas data
continental_sf <- process_data(edge_continental, fuse_continental)
highseas_sf <- process_data(edge_highseas, fuse_highseas)

# Combine continental and high seas data for the combined map
combined_sf <- rbind(continental_sf, highseas_sf)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create color palette
map_pal_raw <- bi_pal(pal = 'PurpleOr', dim = 4, preview = FALSE)
map_pal_mtx <- matrix(map_pal_raw, nrow = 4, ncol = 4)
map_pal_mtx[3, ] <- colorspace::lighten(map_pal_mtx[3, ], .1)
map_pal_mtx[2, ] <- colorspace::lighten(map_pal_mtx[2, ], .2)
map_pal_mtx[1, ] <- colorspace::lighten(map_pal_mtx[1, ], .3)
map_pal_mtx[ , 3] <- colorspace::lighten(map_pal_mtx[ , 3], .1)
map_pal_mtx[ , 2] <- colorspace::lighten(map_pal_mtx[ , 2], .2)
map_pal_mtx[ , 1] <- colorspace::lighten(map_pal_mtx[ , 1], .3)
map_pal_mtx[1, 1] <- '#ffffee'
map_pal <- as.vector(map_pal_mtx) %>% setNames(names(map_pal_raw))

# Color mapping function
get_color <- function(edge, fuse) {
  edge_class <- cut(edge, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  fuse_class <- cut(fuse, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  return(map_pal[(as.numeric(fuse_class)-1)*4 + as.numeric(edge_class)])
}

# Apply colors to all datasets
continental_sf$new_color <- mapply(get_color, continental_sf$EDGE_Priority_Norm, continental_sf$FUSE_Priority_Norm)
highseas_sf$new_color <- mapply(get_color, highseas_sf$EDGE_Priority_Norm, highseas_sf$FUSE_Priority_Norm)
combined_sf$new_color <- mapply(get_color, combined_sf$EDGE_Priority_Norm, combined_sf$FUSE_Priority_Norm)

# Create plot function
create_bivariate_plot <- function(data_sf, title) {
  ggplot() +
    geom_sf(data = data_sf, aes(color = new_color), size = 0.1, alpha = 1) +
    geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
    geom_sf(data = globe_border, fill = NA, color = "grey70", size = 0.5) +
    scale_color_identity() +
    coord_sf() +
    theme_minimal() +
    labs(title = title,
         x = NULL, y = NULL) +
    theme(plot.title = element_text(hjust = 0.5))
}

# Create legend
legend_plot <- bi_legend(pal = map_pal, dim = 4,
                        xlab = 'EDGE2',
                        ylab = 'FUSE')

# Create the individual plots
continental_bivariate <- create_bivariate_plot(continental_sf, "Continental Waters: EDGE2 vs FUSE Priorities")
highseas_bivariate <- create_bivariate_plot(highseas_sf, "High Seas: EDGE2 vs FUSE Priorities")
combined_bivariate_03 <- create_bivariate_plot(combined_sf, "Combined Waters: EDGE2 vs FUSE Priorities; Budget: 0.3")

# Create the map for the ms
# Modify the create_bivariate_plot function for this specific case
combined_bivariate_03_ms <- ggplot() +
    geom_sf(data = combined_sf, aes(color = new_color), size = 0.1, alpha = 1) +
    geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
    geom_sf(data = globe_border, fill = NA, color = "grey70", size = 0.5) +
    scale_color_identity() +
    coord_sf() +
    theme_minimal() +
    labs(x = NULL, y = NULL) +
    annotate("text", x = -Inf, y = Inf, label = "(B)", 
             hjust = -1, vjust = 2, size = 6, fontface = "bold") +
    theme(panel.grid = element_blank(),
          plot.margin = margin(t = 10, r = 10, b = 10, l = 10, unit = "pt"))

# Arrange plots with shared legend
layout <- rbind(
  c(1, 1, 1),
  c(2, 2, 2),
  c(3, 3, 3),
  c(4, 4, 4)
)

# Create legend with larger text
legend_plot <- bi_legend(pal = map_pal, dim = 4,
                        xlab = 'EDGE2',
                        ylab = 'FUSE',
                        size = 2) + # Base size for the legend elements
  theme(
    axis.title = element_text(size = 18, face = "bold"), # Larger axis titles
    axis.text = element_blank(),  # Larger axis text
    legend.text = element_text(size = 12), # Larger legend text
    legend.title = element_text(size = 14, face = "bold") # Larger legend title
  )

# Keep the rest of your grid.arrange code the same
combined_plot <- grid.arrange(
  continental_bivariate,
  highseas_bivariate,
  combined_bivariate_03,
  legend_plot,
  layout_matrix = layout,
  heights = c(0.32, 0.32, 0.32, 0.15),
  widths = unit(c(15, 15, 15), "inches"),
  top = textGrob("Bivariate Maps of EDGE2 and FUSE Priorities", 
                 gp = gpar(fontsize = 16, font = 2))
)

Code
# Display the combined plot
#print(combined_plot)

# Save the plot if needed
# ggsave("bivariate_priority_maps_all.png", combined_plot, width = 15, height = 16, dpi = 300)

Budget: 0.3 with more bins

Code
# Read all RDS files
edge_continental <- readRDS(here::here("Data/EDGE2_full_results_continental_averaged_budget0.3_replicates10.rds"))
edge_highseas <- readRDS(here::here("Data/EDGE2_full_highseas_results_averaged_budget0.3_replicates10.rds"))
fuse_continental <- readRDS(here::here("Data/FUSE_continental_full_results_averaged_budget0.3_replicates10.rds"))
fuse_highseas <- readRDS(here::here("Data/FUSE_full_highseas_results_averaged_budget0.3_replicates10.rds"))

# Get world map data and set projection
world <- ne_countries(scale = "medium", returnclass = "sf")
mcbryde_thomas_2 <- "+proj=mbt_s"

# Function to process and combine data
process_data <- function(edge_data, fuse_data) {
  combined_data <- edge_data %>%
    rename(EDGE_Priority = Priority) %>%
    left_join(fuse_data %>% rename(FUSE_Priority = Priority),
              by = c("Longitude", "Latitude"))
  
  # Normalize priorities to 0-1 range
  combined_data <- combined_data %>%
    mutate(
      EDGE_Priority_Norm = (EDGE_Priority - min(EDGE_Priority)) / (max(EDGE_Priority) - min(EDGE_Priority)),
      FUSE_Priority_Norm = (FUSE_Priority - min(FUSE_Priority)) / (max(FUSE_Priority) - min(FUSE_Priority))
    )
  
  # Transform to sf object
  data_sf <- st_as_sf(combined_data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
    st_transform(crs = mcbryde_thomas_2)
  
  return(data_sf)
}

# Process continental and high seas data
continental_sf <- process_data(edge_continental, fuse_continental)
highseas_sf <- process_data(edge_highseas, fuse_highseas)

# Combine continental and high seas data for the combined map
combined_sf <- rbind(continental_sf, highseas_sf)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create custom 8x8 color palette function
create_custom_bivariate_palette <- function(n = 8) {
  # Create color gradients
  purple_colors <- colorRampPalette(c("#f7f7f7", "#af8dc3", "#7b3294"))(n)
  orange_colors <- colorRampPalette(c("#f7f7f7", "#fdb863", "#e66101"))(n)
  
  # Create the matrix
  pal_matrix <- matrix(NA, nrow = n, ncol = n)
  
  # Fill the matrix with color blends
  for(i in 1:n) {
    for(j in 1:n) {
      # Blend the colors
      color1 <- col2rgb(purple_colors[i])
      color2 <- col2rgb(orange_colors[j])
      
      # Mix the colors with varying weights
      mixed_color <- rgb(
        (color1[1] + color2[1])/2,
        (color1[2] + color2[2])/2,
        (color1[3] + color2[3])/2,
        maxColorValue = 255
      )
      
      pal_matrix[i,j] <- mixed_color
    }
  }
  
  # Apply lightening effect
  for(i in 1:n) {
    for(j in 1:n) {
      light_factor <- (n - i) * 0.1 + (n - j) * 0.1
      if(i != n || j != n) {  # Don't lighten the darkest corner
        pal_matrix[i,j] <- colorspace::lighten(pal_matrix[i,j], light_factor)
      }
    }
  }
  
  # Set the lightest corner
  pal_matrix[1,1] <- "#ffffee"
  
  # Convert matrix to vector
  pal_vector <- as.vector(pal_matrix)
  names(pal_vector) <- paste0("c", 1:(n*n))
  
  return(list(
    palette = pal_vector,
    matrix = pal_matrix
  ))
}

# Create the custom palette
custom_pal <- create_custom_bivariate_palette(8)
map_pal <- custom_pal$palette
map_pal_mtx <- custom_pal$matrix

# Modified color mapping function for 8 bins
get_color <- function(edge, fuse) {
    edge_breaks <- seq(0, 1, length.out = 9)  # Creates 8 bins
    fuse_breaks <- seq(0, 1, length.out = 9)  # Creates 8 bins
    
    edge_class <- cut(edge, breaks = edge_breaks, labels = 1:8, include.lowest = TRUE)
    fuse_class <- cut(fuse, breaks = fuse_breaks, labels = 1:8, include.lowest = TRUE)
    
    return(map_pal[(as.numeric(fuse_class)-1)*8 + as.numeric(edge_class)])
}

# Apply colors to all datasets
continental_sf$new_color <- mapply(get_color, continental_sf$EDGE_Priority_Norm, continental_sf$FUSE_Priority_Norm)
highseas_sf$new_color <- mapply(get_color, highseas_sf$EDGE_Priority_Norm, highseas_sf$FUSE_Priority_Norm)
combined_sf$new_color <- mapply(get_color, combined_sf$EDGE_Priority_Norm, combined_sf$FUSE_Priority_Norm)

# Create plot function
create_bivariate_plot <- function(data_sf, title) {
  ggplot() +
    geom_sf(data = data_sf, aes(color = new_color), size = 0.1, alpha = 1) +
    geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
    geom_sf(data = globe_border, fill = NA, color = "grey70", size = 0.5) +
    scale_color_identity() +
    coord_sf() +
    theme_minimal() +
    labs(title = title,
         x = NULL, y = NULL) +
    theme(plot.title = element_text(hjust = 0.5))
}

# Create the individual plots
continental_bivariate <- create_bivariate_plot(continental_sf, "Continental Waters: EDGE2 vs FUSE Priorities")
highseas_bivariate <- create_bivariate_plot(highseas_sf, "High Seas: EDGE2 vs FUSE Priorities")
combined_bivariate_03 <- create_bivariate_plot(combined_sf, "Combined Waters: EDGE2 vs FUSE Priorities; Budget: 0.3")

# Create the map for the ms
combined_bivariate_03_ms <- ggplot() +
    geom_sf(data = combined_sf, aes(color = new_color), size = 0.1, alpha = 1) +
    geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
    geom_sf(data = globe_border, fill = NA, color = "grey70", size = 0.5) +
    scale_color_identity() +
    coord_sf() +
    theme_minimal() +
    labs(x = NULL, y = NULL) +
    annotate("text", x = -Inf, y = Inf, label = "(B)", 
             hjust = -1, vjust = 2, size = 6, fontface = "bold") +
    theme(panel.grid = element_blank(),
          plot.margin = margin(t = 10, r = 10, b = 10, l = 10, unit = "pt"))

# Create custom legend function for 8x8 grid
create_custom_legend <- function() {
    # Create a data frame for the legend
    legend_data <- expand.grid(
        x = 1:8,
        y = 1:8
    )
    legend_data$color <- as.vector(map_pal_mtx)
    
    # Create the legend plot
    legend_plot <- ggplot(legend_data, aes(x = x, y = y, fill = color)) +
        geom_tile() +
        scale_fill_identity() +
        labs(x = "EDGE2", y = "FUSE") +
        theme_minimal() +
        theme(
            axis.title = element_text(size = 12, face = "bold"),
            axis.text = element_blank(),
            panel.grid = element_blank(),
            plot.margin = margin(t = 5, r = 5, b = 5, l = 5),
            aspect.ratio = 1  # Force square aspect ratio
        ) +
        coord_fixed()  # Ensure square tiles
    
    return(legend_plot)
}

# Create the legend
legend_plot <- create_custom_legend()

# Layout for plots with smaller legend
layout <- rbind(
    c(1, 1, 1, 1),
    c(2, 2, 2, 2),
    c(3, 3, 3, 3),
    c(NA, 4, 4, NA)  # This centers the legend and makes it smaller
)

# Arrange plots with shared legend
combined_plot <- grid.arrange(
    continental_bivariate,
    highseas_bivariate,
    combined_bivariate_03,
    legend_plot,
    layout_matrix = layout,
    heights = c(0.3, 0.3, 0.3, 0.15),
    widths = c(0.25, 0.25, 0.25, 0.25),
    top = textGrob("Bivariate Maps of EDGE2 and FUSE Priorities", 
                   gp = gpar(fontsize = 16, font = 2))
)

Code
# Display the combined plot
print(combined_plot)
TableGrob (5 x 4) "arrange": 5 grobs
  z     cells    name                grob
1 1 (2-2,1-4) arrange      gtable[layout]
2 2 (3-3,1-4) arrange      gtable[layout]
3 3 (4-4,1-4) arrange      gtable[layout]
4 4 (5-5,2-3) arrange      gtable[layout]
5 5 (1-1,1-4) arrange text[GRID.text.665]
Code
# Save the plot if needed
ggsave("bivariate_priority_maps_all_8bins.png", combined_plot, width = 15, height = 16, dpi = 300)

Budget: 0.1

Code
# Read all RDS files
edge_continental <- readRDS(here::here("Data/EDGE2_full_results_continental_averaged_budget0.1_replicates10.rds"))
edge_highseas <- readRDS(here::here("Data/EDGE2_full_highseas_results_averaged_budget0.1_replicates10.rds"))
fuse_continental <- readRDS(here::here("Data/FUSE_full_results_continental_averaged_budget0.1_replicates10.rds"))
fuse_highseas <- readRDS(here::here("Data/FUSE_full_highseas_results_averaged_budget0.1_replicates10.rds"))

# Get world map data and set projection
world <- ne_countries(scale = "medium", returnclass = "sf")
mcbryde_thomas_2 <- "+proj=mbt_s"

# Function to process and combine data
process_data <- function(edge_data, fuse_data) {
  combined_data <- edge_data %>%
    rename(EDGE_Priority = Priority) %>%
    left_join(fuse_data %>% rename(FUSE_Priority = Priority),
              by = c("Longitude", "Latitude"))
  
  # Normalize priorities to 0-1 range
  combined_data <- combined_data %>%
    mutate(
      EDGE_Priority_Norm = (EDGE_Priority - min(EDGE_Priority)) / (max(EDGE_Priority) - min(EDGE_Priority)),
      FUSE_Priority_Norm = (FUSE_Priority - min(FUSE_Priority)) / (max(FUSE_Priority) - min(FUSE_Priority))
    )
  
  # Transform to sf object
  data_sf <- st_as_sf(combined_data, coords = c("Longitude", "Latitude"), crs = 4326) %>%
    st_transform(crs = mcbryde_thomas_2)
  
  return(data_sf)
}

# Process continental and high seas data
continental_sf <- process_data(edge_continental, fuse_continental)
highseas_sf <- process_data(edge_highseas, fuse_highseas)

# Combine continental and high seas data for the combined map
combined_sf <- rbind(continental_sf, highseas_sf)

# Project the world map
world_projected <- st_transform(world, crs = mcbryde_thomas_2)

# Create color palette
map_pal_raw <- bi_pal(pal = 'PurpleOr', dim = 4, preview = FALSE)
map_pal_mtx <- matrix(map_pal_raw, nrow = 4, ncol = 4)
map_pal_mtx[3, ] <- colorspace::lighten(map_pal_mtx[3, ], .1)
map_pal_mtx[2, ] <- colorspace::lighten(map_pal_mtx[2, ], .2)
map_pal_mtx[1, ] <- colorspace::lighten(map_pal_mtx[1, ], .3)
map_pal_mtx[ , 3] <- colorspace::lighten(map_pal_mtx[ , 3], .1)
map_pal_mtx[ , 2] <- colorspace::lighten(map_pal_mtx[ , 2], .2)
map_pal_mtx[ , 1] <- colorspace::lighten(map_pal_mtx[ , 1], .3)
map_pal_mtx[1, 1] <- '#ffffee'
map_pal <- as.vector(map_pal_mtx) %>% setNames(names(map_pal_raw))

# Color mapping function
get_color <- function(edge, fuse) {
  edge_class <- cut(edge, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  fuse_class <- cut(fuse, breaks = c(-Inf, 0.25, 0.5, 0.75, Inf), labels = 1:4)
  return(map_pal[(as.numeric(fuse_class)-1)*4 + as.numeric(edge_class)])
}

# Apply colors to all datasets
continental_sf$new_color <- mapply(get_color, continental_sf$EDGE_Priority_Norm, continental_sf$FUSE_Priority_Norm)
highseas_sf$new_color <- mapply(get_color, highseas_sf$EDGE_Priority_Norm, highseas_sf$FUSE_Priority_Norm)
combined_sf$new_color <- mapply(get_color, combined_sf$EDGE_Priority_Norm, combined_sf$FUSE_Priority_Norm)

# Create plot function
create_bivariate_plot <- function(data_sf, title) {
  ggplot() +
    geom_sf(data = data_sf, aes(color = new_color), size = 0.1, alpha = 1) +
    geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
    geom_sf(data = globe_border, fill = NA, color = "grey70", size = 0.5) +
    scale_color_identity() +
    coord_sf() +
    theme_minimal() +
    labs(title = title,
         x = NULL, y = NULL) +
    theme(plot.title = element_text(hjust = 0.5))
}

# Create legend
legend_plot <- bi_legend(pal = map_pal, dim = 4,
                        xlab = 'EDGE2',
                        ylab = 'FUSE')

# Create the individual plots
continental_bivariate <- create_bivariate_plot(continental_sf, "Continental Waters: EDGE2 vs FUSE Priorities")
highseas_bivariate <- create_bivariate_plot(highseas_sf, "High Seas: EDGE2 vs FUSE Priorities")
combined_bivariate_01 <- create_bivariate_plot(combined_sf, "Combined Waters: EDGE2 vs FUSE Priorities; Budget: 0.1")

#Map for the manuscript 
# Modify the create_bivariate_plot function for this specific case
combined_bivariate_01_ms <- ggplot() +
    geom_sf(data = combined_sf, aes(color = new_color), size = 0.1, alpha = 1) +
    geom_sf(data = world_projected, fill = "lightgray", color = "lightgray") +
    geom_sf(data = globe_border, fill = NA, color = "grey70", size = 0.5) +
    scale_color_identity() +
    coord_sf() +
    theme_minimal() +
    labs(x = NULL, y = NULL) +
    annotate("text", x = -Inf, y = Inf, label = "(A)", 
             hjust = -1, vjust = 2, size = 6, fontface = "bold") +
    theme(panel.grid = element_blank(),
          plot.margin = margin(t = 10, r = 10, b = 10, l = 10, unit = "pt"))

# Arrange plots with shared legend
layout <- rbind(
  c(1, 1, 1),
  c(2, 2, 2),
  c(3, 3, 3),
  c(4, 4, 4)
)

# Create legend with larger text
legend_plot <- bi_legend(pal = map_pal, dim = 4,
                        xlab = 'EDGE2',
                        ylab = 'FUSE',
                        size = 2) + # Base size for the legend elements
  theme(
    axis.title = element_text(size = 18, face = "bold"), # Larger axis titles
    axis.text = element_blank(),  # Larger axis text
    legend.text = element_text(size = 12), # Larger legend text
    legend.title = element_text(size = 14, face = "bold") # Larger legend title
  )

# Keep the rest of your grid.arrange code the same
combined_plot <- grid.arrange(
  continental_bivariate,
  highseas_bivariate,
  combined_bivariate_01,
  legend_plot,
  layout_matrix = layout,
  heights = c(0.32, 0.32, 0.32, 0.15),
  widths = unit(c(15, 15, 15), "inches"),
  top = textGrob("Bivariate Maps of EDGE2 and FUSE Priorities", 
                 gp = gpar(fontsize = 16, font = 2))
)

Code
# Display the combined plot
#print(combined_plot)

# Save the plot if needed
# ggsave("bivariate_priority_maps_all.png", combined_plot, width = 15, height = 16, dpi = 300)

Manuscript maps

Code
layout <- rbind(
  c(1, 1, 1),
  c(2, 2, 2),
  c(3, 3, 3)
)

combined_plot <- grid.arrange(
  combined_bivariate_01_ms,
  combined_bivariate_03_ms,
  legend_plot,
  layout_matrix = layout,
  heights = c(0.32, 0.32, 0.15),
  widths = unit(c(15, 15, 15), "inches") #,
 # top = textGrob("Bivariate Maps of EDGE2 and FUSE Priorities", 
 #                 gp = gpar(fontsize = 16, font = 2))
)

Code
# TIFF version
ggsave(here::here("bivariate_maps_comparison_ms.png"), 
       combined_plot,
       width = 10, 
       height = 12, 
       dpi = 300,
       bg = "white")

#Suplementary figure:
layout <- rbind(
  c(1, 2),
  c(3, 4)
)

combined_plot <- grid.arrange(
  combined_plot_FUSE_01,
  combined_plot_FUSE_03,
  combined_plot_EDGE2_01,
  combined_plot_EDGE2_03,
  layout_matrix = layout,
  top = textGrob("Global Conservation Priorities", 
                 gp = gpar(fontsize = 16, font = 2))
)

Code
# Save if needed
ggsave(here::here("priority_maps_grid.png"), 
       combined_plot,
       width = 15, 
       height = 12, 
       dpi = 300,
       bg = "white")

FUSE : species level priorities

Budget: 0.3

Code
# Protection fraction summary
# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_FUSE_03_continental.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
sp_in_data <- read_csv(here("Data", "continental_puvsp_harmonised.csv"))

# Extract all species names and FUSE values from sp
all_species <- sp$FUSE$info$Species
all_FUSE <- sp$FUSE$info$FUSE

# Get unique species in your data
Species_in_data <- sort(unique(sp_in_data$species_name))

# Create a mapping between species and their FUSE values
species_FUSE_map <- data.frame(
  Species = all_species,
  FUSE = as.numeric(all_FUSE)
)

# Filter the mapping to only include species in your data
filtered_species_FUSE <- species_FUSE_map[species_FUSE_map$Species %in% Species_in_data, ]

# Add Species and FUSE to prot_frac
prot_frac$Species <- filtered_species_FUSE$Species
prot_frac$FUSE <- filtered_species_FUSE$FUSE

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for EDGE2
hist_fuse <- ggplot(prot_frac, aes(x = FUSE)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of FUSE Scores",
       x = "FUSE Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = FUSE, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: FUSE vs Mean Protect Fraction",
       x = "FUSE Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_fuse, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
#High seas waters
# Protection fraction summary for high seas
# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_FUSE_03_highseas.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
sp_in_data <- read_csv(here("Data", "highseas_puvsp_harmonised.csv"))

# Extract all species names and FUSE values from sp
all_species <- sp$FUSE$info$Species
all_FUSE <- sp$FUSE$info$FUSE

# Get unique species in your data
Species_in_data <- sort(unique(sp_in_data$species_name))

# Create a mapping between species and their FUSE values
species_FUSE_map <- data.frame(
  Species = all_species,
  FUSE = as.numeric(all_FUSE)
)

# Filter the mapping to only include species in your data
filtered_species_FUSE <- species_FUSE_map[species_FUSE_map$Species %in% Species_in_data, ]

# Add Species and FUSE to prot_frac
prot_frac$Species <- filtered_species_FUSE$Species
prot_frac$FUSE <- filtered_species_FUSE$FUSE

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction (High Seas)",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for FUSE
hist_fuse <- ggplot(prot_frac, aes(x = FUSE)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of FUSE Scores (High Seas)",
       x = "FUSE Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = FUSE, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: FUSE vs Mean Protect Fraction (High Seas)",
       x = "FUSE Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_fuse, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
#Now combine both and weigth by range size
library(tidyverse)
library(gridExtra)
library(jsonlite)
library(here)

# Load all required data
continental_prot_frac <- readRDS(here::here("Data/protect_fraction_summary_FUSE_03_continental.rds"))
highseas_prot_frac <- readRDS(here::here("Data/protect_fraction_summary_FUSE_03_highseas.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
continental_sp_data <- read_csv(here("Data", "continental_puvsp_harmonised.csv"))
highseas_sp_data <- read_csv(here("Data", "highseas_puvsp_harmonised.csv"))

# Calculate continental range sizes
continental_ranges <- continental_sp_data %>%
  group_by(species_name) %>%
  summarise(continental_range = n())

# Calculate high seas range sizes
highseas_ranges <- highseas_sp_data %>%
  group_by(species_name) %>%
  summarise(highseas_range = n())

# Get species lists
continental_species <- sort(unique(continental_sp_data$species_name))
highseas_species <- sort(unique(highseas_sp_data$species_name))

# Create species-FUSE mapping
all_species <- sp$FUSE$info$Species
all_FUSE <- sp$FUSE$info$FUSE
species_FUSE_map <- data.frame(
  Species = all_species,
  FUSE = as.numeric(all_FUSE)
)

# Add species names to continental data
filtered_continental_FUSE <- species_FUSE_map[species_FUSE_map$Species %in% continental_species, ]
continental_prot_frac$Species <- filtered_continental_FUSE$Species

# Add species names to highseas data
filtered_highseas_FUSE <- species_FUSE_map[species_FUSE_map$Species %in% highseas_species, ]
highseas_prot_frac$Species <- filtered_highseas_FUSE$Species

# Combine the protection fractions with range sizes
combined_protection_FUSE_03 <- full_join(
  continental_prot_frac %>% 
    select(Species, Mean_Protect_Fraction) %>%
    rename(continental_protection = Mean_Protect_Fraction),
  highseas_prot_frac %>% 
    select(Species, Mean_Protect_Fraction) %>%
    rename(highseas_protection = Mean_Protect_Fraction),
  by = "Species"
) %>%
  # Join with the range sizes
  left_join(continental_ranges, by = c("Species" = "species_name")) %>%
  left_join(highseas_ranges, by = c("Species" = "species_name"))

# Calculate weighted protection
combined_protection_FUSE_03 <- combined_protection_FUSE_03 %>%
  mutate(
    # Replace NA with 0 for protection values and ranges
    continental_protection = replace_na(continental_protection, 0),
    highseas_protection = replace_na(highseas_protection, 0),
    continental_range = replace_na(continental_range, 0),
    highseas_range = replace_na(highseas_range, 0),
    # Calculate total range
    total_range = continental_range + highseas_range,
    # Calculate weighted protection
    weighted_protection = (continental_protection * continental_range + 
                         highseas_protection * highseas_range) / 
                         total_range
  )

# Add FUSE scores
combined_protection_FUSE_03 <- left_join(combined_protection_FUSE_03, species_FUSE_map, by = "Species")

# Create summary statistics
summary_stats <- combined_protection_FUSE_03 %>%
  select(-Species) %>%  
  summarise(across(everything(), list(
    min = ~min(., na.rm = TRUE),
    q1 = ~quantile(., 0.25, na.rm = TRUE),
    median = ~median(., na.rm = TRUE),
    mean = ~mean(., na.rm = TRUE),
    q3 = ~quantile(., 0.75, na.rm = TRUE),
    max = ~max(., na.rm = TRUE)
  ))) %>%
  pivot_longer(everything(), 
               names_to = c("variable", "stat"), 
               names_pattern = "(.*)_(.*)") %>%
  pivot_wider(names_from = stat, values_from = value)

# Create and format the flextable
library(flextable)

summary_table <- flextable(summary_stats) %>%
  set_header_labels(
    variable = "Variable",
    min = "Minimum",
    q1 = "1st Quartile",
    median = "Median",
    mean = "Mean",
    q3 = "3rd Quartile",
    max = "Maximum"
  ) %>%
  colformat_double(digits = 3) %>%
  theme_vanilla() %>%
  autofit()

# Display the table
summary_table

Variable

Minimum

1st Quartile

Median

Mean

3rd Quartile

Maximum

continental_protection

0.000

0.302

0.308

0.352

0.336

1.000

highseas_protection

0.000

0.000

0.000

0.101

0.000

1.000

continental_range

0.000

56.000

193.000

1,248.928

650.000

40,875.000

highseas_range

0.000

0.000

0.000

805.935

0.000

63,442.000

total_range

1.000

56.000

200.000

2,054.864

655.000

104,317.000

weighted_protection

0.300

0.303

0.309

0.355

0.341

1.000

FUSE

0.000

0.000

0.001

0.059

0.031

1.000

Code
# Create visualizations
hist_protect <- ggplot(combined_protection_FUSE_03, aes(x = weighted_protection)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Range-Weighted Protection Fraction",
       x = "Weighted Protection Fraction",
       y = "Count")

hist_fuse <- ggplot(combined_protection_FUSE_03, aes(x = FUSE)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of FUSE Scores",
       x = "FUSE Score",
       y = "Count")

scatter_plot <- ggplot(combined_protection_FUSE_03, aes(x = FUSE, y = weighted_protection)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: FUSE vs Weighted Protection Fraction",
       x = "FUSE Score",
       y = "Weighted Protection Fraction")

# Create species range type summary
range_type_summary <- combined_protection_FUSE_03 %>%
  summarise(
    total_species = n(),
    continental_only = sum(highseas_range == 0 & continental_range > 0),
    highseas_only = sum(continental_range == 0 & highseas_range > 0),
    both_ranges = sum(continental_range > 0 & highseas_range > 0)
  ) %>%
  pivot_longer(everything(), 
               names_to = "Distribution Type",
               values_to = "Number of Species") 

# Create and format the flextable
range_type_table <- flextable(range_type_summary) %>%
  set_header_labels(
    `Distribution Type` = "Distribution Type",
    `Number of Species` = "Number of Species"
  ) %>%
  theme_vanilla() %>%
  autofit()

# Display the table
range_type_table

Distribution Type

Number of Species

total_species

1,005

continental_only

802

highseas_only

5

both_ranges

198

Code
# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_fuse, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
# Save the combined protection data
saveRDS(combined_protection_FUSE_03, file = here::here("Data", "combined_protection_FUSE_03.rds"))

Budget: 0.1

Code
# Protection fraction summary
# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_FUSE_01_continental.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
sp_in_data <- read_csv(here("Data", "continental_puvsp_harmonised.csv"))

# Extract all species names and FUSE values from sp
all_species <- sp$FUSE$info$Species
all_FUSE <- sp$FUSE$info$FUSE

# Get unique species in your data
Species_in_data <- sort(unique(sp_in_data$species_name))

# Create a mapping between species and their FUSE values
species_FUSE_map <- data.frame(
  Species = all_species,
  FUSE = as.numeric(all_FUSE)
)

# Filter the mapping to only include species in your data
filtered_species_FUSE <- species_FUSE_map[species_FUSE_map$Species %in% Species_in_data, ]

# Add Species and FUSE to prot_frac
prot_frac$Species <- filtered_species_FUSE$Species
prot_frac$FUSE <- filtered_species_FUSE$FUSE

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction\n(Continental)",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for FUSE
hist_fuse <- ggplot(prot_frac, aes(x = FUSE)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of FUSE Scores\n(Continental)",
       x = "FUSE Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = FUSE, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: FUSE vs Mean Protect Fraction (Continental)",
       x = "FUSE Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_fuse, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
#High seas waters
# Protection fraction summary for high seas
# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_FUSE_01_highseas.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
sp_in_data <- read_csv(here("Data", "highseas_puvsp_harmonised.csv"))

# Extract all species names and FUSE values from sp
all_species <- sp$FUSE$info$Species
all_FUSE <- sp$FUSE$info$FUSE

# Get unique species in your data
Species_in_data <- sort(unique(sp_in_data$species_name))

# Create a mapping between species and their FUSE values
species_FUSE_map <- data.frame(
  Species = all_species,
  FUSE = as.numeric(all_FUSE)
)

# Filter the mapping to only include species in your data
filtered_species_FUSE <- species_FUSE_map[species_FUSE_map$Species %in% Species_in_data, ]

# Add Species and FUSE to prot_frac
prot_frac$Species <- filtered_species_FUSE$Species
prot_frac$FUSE <- filtered_species_FUSE$FUSE

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction\n(High Seas)",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for FUSE
hist_fuse <- ggplot(prot_frac, aes(x = FUSE)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of FUSE Scores\n(High Seas)",
       x = "FUSE Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = FUSE, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: FUSE vs Mean Protect Fraction (High Seas)",
       x = "FUSE Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_fuse, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
#Now combine both and weigth by range size
library(tidyverse)
library(gridExtra)
library(jsonlite)
library(here)

# Load all required data
continental_prot_frac <- readRDS(here::here("Data/protect_fraction_summary_FUSE_01_continental.rds"))
highseas_prot_frac <- readRDS(here::here("Data/protect_fraction_summary_FUSE_01_highseas.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
continental_sp_data <- read_csv(here("Data", "continental_puvsp_harmonised.csv"))
highseas_sp_data <- read_csv(here("Data", "highseas_puvsp_harmonised.csv"))

# Calculate continental range sizes
continental_ranges <- continental_sp_data %>%
  group_by(species_name) %>%
  summarise(continental_range = n())

# Calculate high seas range sizes
highseas_ranges <- highseas_sp_data %>%
  group_by(species_name) %>%
  summarise(highseas_range = n())

# Get species lists
continental_species <- sort(unique(continental_sp_data$species_name))
highseas_species <- sort(unique(highseas_sp_data$species_name))

# Create species-FUSE mapping
all_species <- sp$FUSE$info$Species
all_FUSE <- sp$FUSE$info$FUSE
species_FUSE_map <- data.frame(
  Species = all_species,
  FUSE = as.numeric(all_FUSE)
)

# Add species names to continental data
filtered_continental_FUSE <- species_FUSE_map[species_FUSE_map$Species %in% continental_species, ]
continental_prot_frac$Species <- filtered_continental_FUSE$Species

# Add species names to highseas data
filtered_highseas_FUSE <- species_FUSE_map[species_FUSE_map$Species %in% highseas_species, ]
highseas_prot_frac$Species <- filtered_highseas_FUSE$Species

# Combine the protection fractions with range sizes
combined_protection_FUSE_01 <- full_join(
  continental_prot_frac %>% 
    select(Species, Mean_Protect_Fraction) %>%
    rename(continental_protection = Mean_Protect_Fraction),
  highseas_prot_frac %>% 
    select(Species, Mean_Protect_Fraction) %>%
    rename(highseas_protection = Mean_Protect_Fraction),
  by = "Species"
) %>%
  # Join with the range sizes
  left_join(continental_ranges, by = c("Species" = "species_name")) %>%
  left_join(highseas_ranges, by = c("Species" = "species_name"))

# Calculate weighted protection
combined_protection_FUSE_01 <- combined_protection_FUSE_01 %>%
  mutate(
    # Replace NA with 0 for protection values and ranges
    continental_protection = replace_na(continental_protection, 0),
    highseas_protection = replace_na(highseas_protection, 0),
    continental_range = replace_na(continental_range, 0),
    highseas_range = replace_na(highseas_range, 0),
    # Calculate total range
    total_range = continental_range + highseas_range,
    # Calculate weighted protection
    weighted_protection = (continental_protection * continental_range + 
                         highseas_protection * highseas_range) / 
                         total_range
  )

# Add FUSE scores
combined_protection_FUSE_01 <- left_join(combined_protection_FUSE_01, species_FUSE_map, by = "Species")

# Create summary statistics
summary_stats <- combined_protection_FUSE_01 %>%
  select(-Species) %>%  
  summarise(across(everything(), list(
    min = ~min(., na.rm = TRUE),
    q1 = ~quantile(., 0.25, na.rm = TRUE),
    median = ~median(., na.rm = TRUE),
    mean = ~mean(., na.rm = TRUE),
    q3 = ~quantile(., 0.75, na.rm = TRUE),
    max = ~max(., na.rm = TRUE)
  ))) %>%
  pivot_longer(everything(), 
               names_to = c("variable", "stat"), 
               names_pattern = "(.*)_(.*)") %>%
  pivot_wider(names_from = stat, values_from = value)

# Create and format the flextable
library(flextable)

summary_table <- flextable(summary_stats) %>%
  set_header_labels(
    variable = "Variable",
    min = "Minimum",
    q1 = "1st Quartile",
    median = "Median",
    mean = "Mean",
    q3 = "3rd Quartile",
    max = "Maximum"
  ) %>%
  colformat_double(digits = 3) %>%
  theme_vanilla() %>%
  autofit()

# Display the table
summary_table

Variable

Minimum

1st Quartile

Median

Mean

3rd Quartile

Maximum

continental_protection

0.000

0.103

0.110

0.187

0.194

1.000

highseas_protection

0.000

0.000

0.000

0.075

0.000

1.000

continental_range

0.000

56.000

193.000

1,248.928

650.000

40,875.000

highseas_range

0.000

0.000

0.000

805.935

0.000

63,442.000

total_range

1.000

56.000

200.000

2,054.864

655.000

104,317.000

weighted_protection

0.100

0.103

0.112

0.190

0.200

1.000

FUSE

0.000

0.000

0.001

0.059

0.031

1.000

Code
# Create visualizations
hist_protect <- ggplot(combined_protection_FUSE_01, aes(x = weighted_protection)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Range-Weighted Protection Fraction",
       x = "Weighted Protection Fraction",
       y = "Count")

hist_fuse <- ggplot(combined_protection_FUSE_01, aes(x = FUSE)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of FUSE Scores",
       x = "FUSE Score",
       y = "Count")

scatter_plot <- ggplot(combined_protection_FUSE_01, aes(x = FUSE, y = weighted_protection)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: FUSE vs Weighted Protection Fraction",
       x = "FUSE Score",
       y = "Weighted Protection Fraction")

# Create species range type summary
range_type_summary <- combined_protection_FUSE_01 %>%
  summarise(
    total_species = n(),
    continental_only = sum(highseas_range == 0 & continental_range > 0),
    highseas_only = sum(continental_range == 0 & highseas_range > 0),
    both_ranges = sum(continental_range > 0 & highseas_range > 0)
  ) %>%
  pivot_longer(everything(), 
               names_to = "Distribution Type",
               values_to = "Number of Species") 

# Create and format the flextable
range_type_table <- flextable(range_type_summary) %>%
  set_header_labels(
    `Distribution Type` = "Distribution Type",
    `Number of Species` = "Number of Species"
  ) %>%
  theme_vanilla() %>%
  autofit()

# Display the table
range_type_table

Distribution Type

Number of Species

total_species

1,005

continental_only

802

highseas_only

5

both_ranges

198

Code
# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_fuse, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
# Save the combined protection data
saveRDS(combined_protection_FUSE_01, file = here::here("Data", "combined_protection_FUSE_01.rds"))

EDGE2 : species level priorities

Budget: 0.3

Code
# Protection fraction summary
# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE2_03_continental.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
sp_in_data <- read_csv(here("Data", "continental_puvsp_harmonised.csv"))

# Extract all species names and EDGE2 values from sp
all_species <- sp$EDGE2$info$Species
all_EDGE2 <- sp$EDGE2$info$EDGE2

# Get unique species in your data
Species_in_data <- sort(unique(sp_in_data$species_name))

# Create a mapping between species and their EDGE2 values
species_EDGE2_map <- data.frame(
  Species = all_species,
  EDGE2 = as.numeric(all_EDGE2)
)

# Filter the mapping to only include species in your data
filtered_species_EDGE2 <- species_EDGE2_map[species_EDGE2_map$Species %in% Species_in_data, ]

# Add Species and EDGE2 to prot_frac
prot_frac$Species <- filtered_species_EDGE2$Species
prot_frac$EDGE2 <- filtered_species_EDGE2$EDGE2

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction\n(Continental)",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for EDGE2
hist_EDGE2 <- ggplot(prot_frac, aes(x = EDGE2)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of EDGE2 Scores\n(Continental)",
       x = "EDGE2 Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = EDGE2, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: EDGE2 vs Mean Protect Fraction (Continental)",
       x = "EDGE2 Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_EDGE2, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
#High seas waters
# Protection fraction summary for high seas
# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE2_03_highseas.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
sp_in_data <- read_csv(here("Data", "highseas_puvsp_harmonised.csv"))

# Extract all species names and EDGE2 values from sp
all_species <- sp$EDGE2$info$Species
all_EDGE2 <- sp$EDGE2$info$EDGE2

# Get unique species in your data
Species_in_data <- sort(unique(sp_in_data$species_name))

# Create a mapping between species and their EDGE2 values
species_EDGE2_map <- data.frame(
  Species = all_species,
  EDGE2 = as.numeric(all_EDGE2)
)

# Filter the mapping to only include species in your data
filtered_species_EDGE2 <- species_EDGE2_map[species_EDGE2_map$Species %in% Species_in_data, ]

# Add Species and EDGE2 to prot_frac
prot_frac$Species <- filtered_species_EDGE2$Species
prot_frac$EDGE2 <- filtered_species_EDGE2$EDGE2

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction\n(High Seas)",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for EDGE2
hist_EDGE2 <- ggplot(prot_frac, aes(x = EDGE2)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of EDGE2 Scores\n(High Seas)",
       x = "EDGE2 Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = EDGE2, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: EDGE2 vs Mean Protect Fraction (High Seas)",
       x = "EDGE2 Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_EDGE2, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
#Now combine both and weigth by range size
library(tidyverse)
library(gridExtra)
library(jsonlite)
library(here)

# Load all required data
continental_prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE2_03_continental.rds"))
highseas_prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE2_03_highseas.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
continental_sp_data <- read_csv(here("Data", "continental_puvsp_harmonised.csv"))
highseas_sp_data <- read_csv(here("Data", "highseas_puvsp_harmonised.csv"))

# Calculate continental range sizes
continental_ranges <- continental_sp_data %>%
  group_by(species_name) %>%
  summarise(continental_range = n())

# Calculate high seas range sizes
highseas_ranges <- highseas_sp_data %>%
  group_by(species_name) %>%
  summarise(highseas_range = n())

# Get species lists
continental_species <- sort(unique(continental_sp_data$species_name))
highseas_species <- sort(unique(highseas_sp_data$species_name))

# Create species-EDGE2 mapping
all_species <- sp$EDGE2$info$Species
all_EDGE2 <- sp$EDGE2$info$EDGE2
species_EDGE2_map <- data.frame(
  Species = all_species,
  EDGE2 = as.numeric(all_EDGE2)
)

# Add species names to continental data
filtered_continental_EDGE2 <- species_EDGE2_map[species_EDGE2_map$Species %in% continental_species, ]
continental_prot_frac$Species <- filtered_continental_EDGE2$Species

# Add species names to highseas data
filtered_highseas_EDGE2 <- species_EDGE2_map[species_EDGE2_map$Species %in% highseas_species, ]
highseas_prot_frac$Species <- filtered_highseas_EDGE2$Species

# Combine the protection fractions with range sizes
combined_protection_EDGE2_03 <- full_join(
  continental_prot_frac %>% 
    select(Species, Mean_Protect_Fraction) %>%
    rename(continental_protection = Mean_Protect_Fraction),
  highseas_prot_frac %>% 
    select(Species, Mean_Protect_Fraction) %>%
    rename(highseas_protection = Mean_Protect_Fraction),
  by = "Species"
) %>%
  # Join with the range sizes
  left_join(continental_ranges, by = c("Species" = "species_name")) %>%
  left_join(highseas_ranges, by = c("Species" = "species_name"))

# Calculate weighted protection
combined_protection_EDGE2_03 <- combined_protection_EDGE2_03 %>%
  mutate(
    # Replace NA with 0 for protection values and ranges
    continental_protection = replace_na(continental_protection, 0),
    highseas_protection = replace_na(highseas_protection, 0),
    continental_range = replace_na(continental_range, 0),
    highseas_range = replace_na(highseas_range, 0),
    # Calculate total range
    total_range = continental_range + highseas_range,
    # Calculate weighted protection
    weighted_protection = (continental_protection * continental_range + 
                         highseas_protection * highseas_range) / 
                         total_range
  )

# Add EDGE2 scores
combined_protection_EDGE2_03 <- left_join(combined_protection_EDGE2_03, species_EDGE2_map, by = "Species")

# Create summary statistics
summary_stats <- combined_protection_EDGE2_03 %>%
  select(-Species) %>%  
  summarise(across(everything(), list(
    min = ~min(., na.rm = TRUE),
    q1 = ~quantile(., 0.25, na.rm = TRUE),
    median = ~median(., na.rm = TRUE),
    mean = ~mean(., na.rm = TRUE),
    q3 = ~quantile(., 0.75, na.rm = TRUE),
    max = ~max(., na.rm = TRUE)
  ))) %>%
  pivot_longer(everything(), 
               names_to = c("variable", "stat"), 
               names_pattern = "(.*)_(.*)") %>%
  pivot_wider(names_from = stat, values_from = value)

# Create and format the flextable
library(flextable)

summary_table <- flextable(summary_stats) %>%
  set_header_labels(
    variable = "Variable",
    min = "Minimum",
    q1 = "1st Quartile",
    median = "Median",
    mean = "Mean",
    q3 = "3rd Quartile",
    max = "Maximum"
  ) %>%
  colformat_double(digits = 3) %>%
  theme_vanilla() %>%
  autofit()

# Display the table
summary_table

Variable

Minimum

1st Quartile

Median

Mean

3rd Quartile

Maximum

continental_protection

0.000

0.302

0.307

0.348

0.330

1.000

highseas_protection

0.000

0.000

0.000

0.096

0.000

1.000

continental_range

0.000

56.000

193.000

1,248.928

650.000

40,875.000

highseas_range

0.000

0.000

0.000

805.935

0.000

63,442.000

total_range

1.000

56.000

200.000

2,054.864

655.000

104,317.000

weighted_protection

0.300

0.303

0.308

0.351

0.333

1.000

EDGE2

0.000

0.000

0.001

0.045

0.018

1.000

Code
# Create visualizations
hist_protect <- ggplot(combined_protection_EDGE2_03, aes(x = weighted_protection)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Range-Weighted Protection Fraction",
       x = "Weighted Protection Fraction",
       y = "Count")

hist_EDGE2 <- ggplot(combined_protection_EDGE2_03, aes(x = EDGE2)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of EDGE2 Scores",
       x = "EDGE2 Score",
       y = "Count")

scatter_plot <- ggplot(combined_protection_EDGE2_03, aes(x = EDGE2, y = weighted_protection)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: EDGE2 vs Weighted Protection Fraction",
       x = "EDGE2 Score",
       y = "Weighted Protection Fraction")

# Create species range type summary
range_type_summary <- combined_protection_EDGE2_03 %>%
  summarise(
    total_species = n(),
    continental_only = sum(highseas_range == 0 & continental_range > 0),
    highseas_only = sum(continental_range == 0 & highseas_range > 0),
    both_ranges = sum(continental_range > 0 & highseas_range > 0)
  ) %>%
  pivot_longer(everything(), 
               names_to = "Distribution Type",
               values_to = "Number of Species") 

# Create and format the flextable
range_type_table <- flextable(range_type_summary) %>%
  set_header_labels(
    `Distribution Type` = "Distribution Type",
    `Number of Species` = "Number of Species"
  ) %>%
  theme_vanilla() %>%
  autofit()

# Display the table
range_type_table

Distribution Type

Number of Species

total_species

1,005

continental_only

802

highseas_only

5

both_ranges

198

Code
# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_EDGE2, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
# Save the combined protection data
saveRDS(combined_protection_EDGE2_03, file = here::here("Data", "combined_protection_EDGE2_03.rds"))

Budget: 0.1

Code
# Protection fraction summary
# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE2_01_continental.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
sp_in_data <- read_csv(here("Data", "continental_puvsp_harmonised.csv"))

# Extract all species names and EDGE2 values from sp
all_species <- sp$EDGE2$info$Species
all_EDGE2 <- sp$EDGE2$info$EDGE2

# Get unique species in your data
Species_in_data <- sort(unique(sp_in_data$species_name))

# Create a mapping between species and their EDGE2 values
species_EDGE2_map <- data.frame(
  Species = all_species,
  EDGE2 = as.numeric(all_EDGE2)
)

# Filter the mapping to only include species in your data
filtered_species_EDGE2 <- species_EDGE2_map[species_EDGE2_map$Species %in% Species_in_data, ]

# Add Species and EDGE2 to prot_frac
prot_frac$Species <- filtered_species_EDGE2$Species
prot_frac$EDGE2 <- filtered_species_EDGE2$EDGE2

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction\n(Continental)",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for EDGE2
hist_EDGE2 <- ggplot(prot_frac, aes(x = EDGE2)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of EDGE2 Scores\n(Continental)",
       x = "EDGE2 Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = EDGE2, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: EDGE2 vs Mean Protect Fraction (Continental)",
       x = "EDGE2 Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_EDGE2, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
#High seas waters
# Protection fraction summary for high seas
# Read the data
prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE2_01_highseas.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
sp_in_data <- read_csv(here("Data", "highseas_puvsp_harmonised.csv"))

# Extract all species names and EDGE2 values from sp
all_species <- sp$EDGE2$info$Species
all_EDGE2 <- sp$EDGE2$info$EDGE2

# Get unique species in your data
Species_in_data <- sort(unique(sp_in_data$species_name))

# Create a mapping between species and their EDGE2 values
species_EDGE2_map <- data.frame(
  Species = all_species,
  EDGE2 = as.numeric(all_EDGE2)
)

# Filter the mapping to only include species in your data
filtered_species_EDGE2 <- species_EDGE2_map[species_EDGE2_map$Species %in% Species_in_data, ]

# Add Species and EDGE2 to prot_frac
prot_frac$Species <- filtered_species_EDGE2$Species
prot_frac$EDGE2 <- filtered_species_EDGE2$EDGE2

# Create histogram for Mean_Protect_Fraction
hist_protect <- ggplot(prot_frac, aes(x = Mean_Protect_Fraction)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Mean Protect Fraction\n(High Seas)",
       x = "Mean Protect Fraction",
       y = "Count")

# Create histogram for EDGE2
hist_EDGE2 <- ggplot(prot_frac, aes(x = EDGE2)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of EDGE2 Scores\n(High Seas)",
       x = "EDGE2 Score",
       y = "Count")

# Create scatterplot
scatter_plot <- ggplot(prot_frac, aes(x = EDGE2, y = Mean_Protect_Fraction)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: EDGE2 vs Mean Protect Fraction (High Seas)",
       x = "EDGE2 Score",
       y = "Mean Protect Fraction")

# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_EDGE2, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
library(tidyverse)
library(gridExtra)
library(jsonlite)
library(here)

# Load all required data
continental_prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE2_01_continental.rds"))
highseas_prot_frac <- readRDS(here::here("Data/protect_fraction_summary_EDGE2_01_highseas.rds"))
sp <- fromJSON(here("Data", "shark_conservation_metrics_no_freshwater.json"))
continental_sp_data <- read_csv(here("Data", "continental_puvsp_harmonised.csv"))
highseas_sp_data <- read_csv(here("Data", "highseas_puvsp_harmonised.csv"))

# Calculate continental range sizes
continental_ranges <- continental_sp_data %>%
  group_by(species_name) %>%
  summarise(continental_range = n())

# Calculate high seas range sizes
highseas_ranges <- highseas_sp_data %>%
  group_by(species_name) %>%
  summarise(highseas_range = n())

# Get species lists
continental_species <- sort(unique(continental_sp_data$species_name))
highseas_species <- sort(unique(highseas_sp_data$species_name))

# Create species-EDGE2 mapping
all_species <- sp$EDGE2$info$Species
all_EDGE2 <- sp$EDGE2$info$EDGE2
species_EDGE2_map <- data.frame(
  Species = all_species,
  EDGE2 = as.numeric(all_EDGE2)
)

# Add species names to continental data
filtered_continental_EDGE2 <- species_EDGE2_map[species_EDGE2_map$Species %in% continental_species, ]
continental_prot_frac$Species <- filtered_continental_EDGE2$Species

# Add species names to highseas data
filtered_highseas_EDGE2 <- species_EDGE2_map[species_EDGE2_map$Species %in% highseas_species, ]
highseas_prot_frac$Species <- filtered_highseas_EDGE2$Species

# Combine the protection fractions with range sizes
combined_protection_EDGE2_01 <- full_join(
  continental_prot_frac %>% 
    select(Species, Mean_Protect_Fraction) %>%
    rename(continental_protection = Mean_Protect_Fraction),
  highseas_prot_frac %>% 
    select(Species, Mean_Protect_Fraction) %>%
    rename(highseas_protection = Mean_Protect_Fraction),
  by = "Species"
) %>%
  # Join with the range sizes
  left_join(continental_ranges, by = c("Species" = "species_name")) %>%
  left_join(highseas_ranges, by = c("Species" = "species_name"))

# Calculate weighted protection
combined_protection_EDGE2_01 <- combined_protection_EDGE2_01 %>%
  mutate(
    # Replace NA with 0 for protection values and ranges
    continental_protection = replace_na(continental_protection, 0),
    highseas_protection = replace_na(highseas_protection, 0),
    continental_range = replace_na(continental_range, 0),
    highseas_range = replace_na(highseas_range, 0),
    # Calculate total range
    total_range = continental_range + highseas_range,
    # Calculate weighted protection
    weighted_protection = (continental_protection * continental_range + 
                         highseas_protection * highseas_range) / 
                         total_range
  )

# Add EDGE2 scores
combined_protection_EDGE2_01 <- left_join(combined_protection_EDGE2_01, species_EDGE2_map, by = "Species")

# Create summary statistics
summary_stats <- combined_protection_EDGE2_01 %>%
  select(-Species) %>%  
  summarise(across(everything(), list(
    min = ~min(., na.rm = TRUE),
    q1 = ~quantile(., 0.25, na.rm = TRUE),
    median = ~median(., na.rm = TRUE),
    mean = ~mean(., na.rm = TRUE),
    q3 = ~quantile(., 0.75, na.rm = TRUE),
    max = ~max(., na.rm = TRUE)
  ))) %>%
  pivot_longer(everything(), 
               names_to = c("variable", "stat"), 
               names_pattern = "(.*)_(.*)") %>%
  pivot_wider(names_from = stat, values_from = value)

# Create and format the flextable
library(flextable)

summary_table <- flextable(summary_stats) %>%
  set_header_labels(
    variable = "Variable",
    min = "Minimum",
    q1 = "1st Quartile",
    median = "Median",
    mean = "Mean",
    q3 = "3rd Quartile",
    max = "Maximum"
  ) %>%
  colformat_double(digits = 3) %>%
  theme_vanilla() %>%
  autofit()

# Display the table
summary_table

Variable

Minimum

1st Quartile

Median

Mean

3rd Quartile

Maximum

continental_protection

0.000

0.103

0.109

0.174

0.160

1.000

highseas_protection

0.000

0.000

0.000

0.068

0.000

1.000

continental_range

0.000

56.000

193.000

1,248.928

650.000

40,875.000

highseas_range

0.000

0.000

0.000

805.935

0.000

63,442.000

total_range

1.000

56.000

200.000

2,054.864

655.000

104,317.000

weighted_protection

0.100

0.103

0.111

0.177

0.168

1.000

EDGE2

0.000

0.000

0.001

0.045

0.018

1.000

Code
# Create visualizations
hist_protect <- ggplot(combined_protection_EDGE2_01, aes(x = weighted_protection)) +
  geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
  scale_x_continuous(limits=c(0,1)) + 
  theme_minimal() +
  labs(title = "Histogram of Range-Weighted Protection Fraction",
       x = "Weighted Protection Fraction",
       y = "Count")

hist_EDGE2 <- ggplot(combined_protection_EDGE2_01, aes(x = EDGE2)) +
  geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of EDGE2 Scores",
       x = "EDGE2 Score",
       y = "Count")

scatter_plot <- ggplot(combined_protection_EDGE2_01, aes(x = EDGE2, y = weighted_protection)) +
  geom_point(alpha = 0.6, color = "darkblue") +
  theme_minimal() +
  scale_y_continuous(limits=c(0,1)) + 
  labs(title = "Scatterplot: EDGE2 vs Weighted Protection Fraction",
       x = "EDGE2 Score",
       y = "Weighted Protection Fraction")

# Create species range type summary
range_type_summary <- combined_protection_EDGE2_01 %>%
  summarise(
    total_species = n(),
    continental_only = sum(highseas_range == 0 & continental_range > 0),
    highseas_only = sum(continental_range == 0 & highseas_range > 0),
    both_ranges = sum(continental_range > 0 & highseas_range > 0)
  ) %>%
  pivot_longer(everything(), 
               names_to = "Distribution Type",
               values_to = "Number of Species") 

# Create and format the flextable
range_type_table <- flextable(range_type_summary) %>%
  set_header_labels(
    `Distribution Type` = "Distribution Type",
    `Number of Species` = "Number of Species"
  ) %>%
  theme_vanilla() %>%
  autofit()

# Display the table
range_type_table

Distribution Type

Number of Species

total_species

1,005

continental_only

802

highseas_only

5

both_ranges

198

Code
# Arrange plots in a grid
grid_plot <- grid.arrange(
  hist_protect, hist_EDGE2, scatter_plot,
  layout_matrix = rbind(c(1,2), c(3,3)),
  widths = c(1, 1),
  heights = c(1, 1)
)

Code
# Save the combined protection data
saveRDS(combined_protection_EDGE2_01, file = here::here("Data", "combined_protection_EDGE2_01.rds"))